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Abstract. A Green's function based solver for the modified Bessel equation has been developed 
with the primary motivation of solving the Poisson equation in cylindrical geometries. The method 
is implemented using a Discrete Hankel Transform and a Green's function based on the modified 
Bessel functions of the first and second kind. The computation of these Bessel functions has been 
implemented to avoid scaling problems due to their exponential and singular behavior, allowing the 
method to be used for large order problems, as would arise in solving the Poisson equation with a 
dense azimuthal grid. The method has been tested on monotonically decaying and oscillatory inputs, 
checking for errors due to interpolation and/or aliasing. The error has been found to reach machine 
precision and to have computational time linearly proportional to the number of nodes. 
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1. Introduction. This paper is motivated by the requirement for a Poisson 
solver for cyhndrical domains. Such a solver is a basic element in solving a range of 
physical problems and accurate methods have been developed for Cartesian [5] and 
cylindrical [3] domains. An issue which arises in solving the problem in a cylindrical 
coordinate system is the singularity which arises at the axis due to the form of the 
differential operator. A recent paper by Pataki and Greengard [11] introduces a 
Green's function solver for the Poisson equation, which avoids problems with this 
singularity, and automatically imposes a radiation boundary condition, by using an 
integral formulation for the solution. 

The method which Pataki and Greengard [11] develop is based on Fourier trans- 
forms in the axial and azimuthal coordinate, followed by the solution of a modified 
Bessel equation in the radial coordinate. They test the accuracy of their method on 
an axisymmetric problem with monotonic decay in radius, and show its application 
to an asymmetric problem. A problem which arises in their algorithm is that the 
integration technique used does not work well for large azimuthal orders, i.e. meshes 
dense in angle, or for large axial wavenumbers, i.e. meshes dense in the axial coor- 
dinate, due to the poor scaling of the modified Bessel functions which appear in the 
Green's function for the problem. 

This motivated an attempt to find a robust method for solving the modified 
Bessel equation, which would work for large wavenumbers and azimuthal orders. The 
method developed, which is described in the rest of this paper, is based on the Discrete 
Hankel Transform (DHT), tabulated integrals, and recursions for the modified Bessel 
functions. 

2. Problem formulation. The Poisson equation in cylindrical coordinates is: 

Urr{r, 0, z) + -Ur{r, 6*, z) + -ugeir, 6*, z) + u^^(r, 9, z) = f{r, 0, z), (2.1) 
r r 

where (r, 9, z) are cylindrical coordinates, u is the solution and / is some forcing term. 
With the problem defined on nodes regularly spaced in 9 and z, this equation can be 
solved [11] by using the FFT to Fourier transform u and f m 9 and z to yield a set 
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of modified Bessel equations: 

+ - +'^') = (2.2) 

where n is the azimuthal order and n the axial wavenumber, with u^") and /^"^ the 
Fourier transformed solution and forcing term respectively. After solving Equation 2.2 
for each value of n and k, (r, k) can be inverse Fourier transformed to yield the 
solution u{r, 9, z). 

Pataki and Grccngard [11] give a method for solving this modified Bcsscl equation, 
subject to a radiation boundary condition at some outer radius r = R. This is done 
using the Green's function for the modified Bessel equation with the solution written: 



The Green's function G„ is: 



/ G„(K,r, s)/(s,K)ds. (2.3) 
Jo 



where /„(a;) and if„(x) are the modified Bessel functions of the first and second kind 
respectively. The first of these, 7„, grows exponentially, while decays exponen- 
tially, but has log x behavior at the origin, leading to some numerical problems 
caused by the scaling of the Green's function. In practice, the solution is computed 
as: 

J\n\KKm) Jq 

- In{i^T)K^{KRm) j" ^TIT^ ds, (2.5) 

where the radial domain is divided at radii Rm which are used to set reference val- 
ues of the modified Bessel functions. The product In{Ks)Kn{KRm), and the ratio 
Kn{ns) / Kn{KRm)i are thus well-scaled avoiding problems in computation, as long as 
|k(s — i?m)| is not too large. Pataki [10] reports that the integration is performed using 
a dyadic grid, and that by scaling the modified Bessel functions as they arc computed, 
the method works well for n < 40, corresponding to 80 points in the azimuthal mesh. 

For many applications, it is desirable to use a denser mesh than this and so a 
different approach was sought for the solution of Equation 2.2. The natural transform 
technique for problems in polar coordinates is the Discrete Hankel Transform (DHT), 
which expands a function as a series of ordinary Bessel functions Jm {x) . For a function 
/(r), 0<r<R: 

M 

fir) « fmJmiamr) (2.6) 

m=l 

where Jm{c(mR) = and denotes the m coefficient of the expansion of /. If 
the DHT of f^"'\r,K) is available, the solution of Equation 2.2 can be immediately 
written: 

u(")(r,K) = V / Jm{ams)G„{K,r,s)ds. (2.7) 

m=l 
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The integrals required in Equation 2.7 are given in standard tables [7]: 

/ In{Ks)Jn{as)sds = [a J„+i (ar)/„ (wr) + k7„+i (Kr)J„(ar)] ^— — 

Jo " + 

+ arJn+i{ar)Kn{nr) - KrKn+i{Kr)Jn{ar) 



Kn{Ks)Jn{cts)s ds 



(2.8a) 

1 

(2.8b) 



which gives a solution for the problem in terms of the DHT coefficients and 
K. As written, this solution is correct, but not numerically useful, due to the poor 
scaling of the modified Bessel functions, especially for large values of k and/or n. It 
must be rewritten in order to avoid numerical difficulties. 

3. Numerical implementation. In order to avoid numerical difficulties caused 
by poor scaling of the modified Bessel function, Equations 2.8 are rewritten and used 
to give the convolution of the Green's function with the ordinary Bessel function as: 



Jo 



Gn{K,r, s)Jn{as) ds 



R 



In{Kr)Kn{KR) 



aJn+i(aR) - KJn{aR) 

Kn[KR) 



■ nr 



In{Kr)Kn{Kr)Jn{ar) 



1 



, r = 0, n — 0, 



I„{Kr) Kninr) 



: 0, r = Q,n^O. 



(3.1a) 

(3.1b) 
(3.1c) 



Written in this form, the modified Bessel functions appear only as ratios /„+i (x) / In (x) 
and Kn+i{x) / Kn{x) or as the products In{x)Kn{x) and In{Kr)Kn(KR). The ratios 
can be computed directly using standard functional relations, while the products are 
calculated using the same ratios combined with modified Bessel fimctions of order 
zero, which can be computed accurately and stably. Implementation of the solution 
technique thus requires two elements, a method for the calculation of ratios of modified 
Bessel functions, and a method for computing the DHT. In practice, the input will not 
be defined on the nodes of the DHT, so an interpolation scheme will also be required. 
The method has been coded making use of the GNU Scientific Library (GSL) [4], 
which provides functions for the computation of scaled versions of the modified Bessel 
functions, directly returning /„(a;) exp[— a;] and Kn{x) exp[a;]. The algorithm has been 
designed to use these scaled functions, to avoid problems of underfiow and overflow. 

3.1. Ratios of modified Bessel functions. The ratios of modified Bessel func- 
tions can be computed using standard functional relations [7, 8.486]: 



2n 

In-l{x) = —In{x)+In+l{x), 

2n 

Kn+l{x) = —Kn{x) + Kn-l{x), 



(3.2a) 

(3.2b) 
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using an approach similar to that of Amos [1], who writes the ratios of successive 
functions as: 



In-l{x) 
Kn+l{x) 
Kn{x) 



X 



2n- 
2n 



+ 



Xln+lix)/ln{x)'' 
Kn-l{x) 
Kn{x) • 



(3.3a) 
(3.3b) 



The recursion for /„ [x) / (cc) is stable for descending n while that for if„+i {x)/Kn{x) 
is stable for increasing n. The recursion for Kn+\ [x) / Kn [x) is seeded with Ki [x) /Kq {x) , 
computed using the scaled form of Ko{x) and Ki{x). The recursion for /„(.x)//„_i(a;) 
is seeded using Olver's asymptotic formula [9, 10.41.10] for modified Bessel functions 
of large order, as recommended by Amos [1] , starting at order equal to the larger of 
n + 8 and 32. The asymptotic expansion is given by: 



E 



(27rn)V2(i + ^2)1/4 ^ nq ' 
x/n,7) = {l + z^f/'^ + \og 



(3.4) 



1+ (1+^2)1/2 



t =1/(1 + ^2)- 



-1/2 



Uo{t) 
Ul{t) 
U2{t) 



= 1, 



(3t-5i^)/24, 

{81t^ - 462t^ + 385i^)/1152, 

(30375i^ - 369603t^ + 765765i^ - 425425t^)/414720, 



while for small arguments, (a;/2)^ < n + 1, the series expansion of In{x) is used [7, 
8.445]. 

Given a sequence of ratios of modified Bessel functions, the products in Equa- 
tion 3.1 can be computed as: 



n-l 



i=0 



li+l 



and 



In{Kr)Kn{KR) : 
A: 



[/o(«:r)e-' 

K{r-R)/n 



[/fo(«i?)e"^] n ^ 



Ii{K,r) 



Ii+i{Kr) 
Ii{Kr) 



Ki+i{Kr) 



Ki{Kr) 



Ki{KR) 



(3.5) 



(3.6) 



where terms of the form /„(.t) cxp[— x] and Kn{x) exp[a;] are computed directly using 
the scaled form of the modified Bessel functions. The ratios of successive modified 
Bessel functions are well scaled and multiplying them in pairs as in the products of 
Equations 3.5 and 3.6 avoids underflow and overflow problems. 

3.2. Discrete Hankel Transform. The coefficients of the DHT arc computed 
using the method of Lcmoinc [8] . This is essentially a quadrature rule based on the 
zeros of the ordinary Bessel function of order n, Jn{x). The function to be transformed 
is specified at these zeros Xm, < m < M, Jn{xm) = 0, and the DHT is given by a 
matrix multiplication of the vector of input data f{xm) with the matrix entries given 
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by: 



Xm \Jn+l{Xm)Jn+l{Xj)\' 



In the calculations presented here, the GSL implementation [4] of Lemoine's method 
was used, but with a modification to compute the zeros of J„(a;) using the 0(M) 
algorithm of Glaser et. al [6]. 

In order to compute the DHT, the input must be specified at the zeros of the 
Bessel function. If only one order n is of interest, this presents no difiiculties, but if 
the solution is to be found for multiple values of n, as in solving a Poisson equation, for 
example, an interpolation scheme is required to transfer the input from the problem 
mesh onto the DHT nodes, in particular because the zeros are not the same for 
different orders of Bessel function. 

3.3. Interpolation. Given that an interpolation scheme will almost always be 
needed, the approach used by Pataki and Greengard [11] has been adopted. The 
domain < r < i? is divided into N blocks, i?„ < r < i?„+i, n = 0, . . . , TV — 1. Each 
block is discretized with P points, given by the Chebyshev nodes of the second kind: 

Rn+l + Rn , Rn+1 — Rn „ , □ /n o\ 

rnP+p = ^ + ^ cos—, p = 0,l,...,P. (3.8) 

Evaluation of /'"^ (r, k) within each block is performed using barycentric Lagrangian 
interpolation [2]. Since Equation 3.1 can be computed directly at arbitrary values of 
r, the solution is generated on the input nodes, with no requirement for interpolation 
from the DHT nodes. 

3.4. Summary of algorithm. Given the elements described above, the solution 
algorithm can be summarized as follows. For a given order n, wavenumber k and input 
f'^"\r,K), 0<r<R: 

1. generate, if necessary, the DHT matrix and corresponding nodes r^; 

2. if necessary, interpolate f^^^r, k) onto the nodes r^; 

3. perform the DHT to yield 

4. evaluate Equation 2.7 at the input nodes using Equations 3.1. 

4. Numerical tests. The solution method is tested using a function which can 
be varied to examine the performance of the algorithm with regard to potential sources 
of error. The main sources of error in the algorithm arise from the interpolation 
schemes and aliasing. These errors arise in both the direct method, where the input 
is specified on the DHT nodes, and when the input must be interpolated from another 
mesh onto these points. 

Interpolation errors arise when the interpolation scheme is unable to accurately 
resolve the function which is being interpolated. This can occur because the interpo- 
lation method proper does not have the required properties to give a well-converged 
estimate of the underlying function, or because the interpolation nodes are not dense 
enough to take advantage of an otherwise good interpolation method. In this sense, 
'interpolation scheme' refers both to the explicitly stated polynomial interpolation 
method used to transfer data from the input mesh to the DHT nodes, and to the 
interpolation which is performed implicitly in the quadrature scheme of the DHT. 

The second source of error is aliasing, when the point distribution is not dense 
enough to capture the spatial frequencies present in the input. This can happen in 
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Fig. 4.1. Test function with a = l, fS = S, m = 8 

the DHT, if the analytically defined input has wavenumbers am with m > M, so that 
the expansion of Equation 2.6 does not contain the full set of radial wavenumbers 
am present in the input. Clearly, it can also happen in the Chebyshev interpolation 
scheme if the node density is not high enough, even if the DHT would otherwise 
contain enough nodes to capture the full behavior of the input. 

In order to assess the performance of the algorithm against these criteria, the test 
function: 

u(")(r) = (r/r„,ax)"e-('^'-''-«'/"' cos/3r, m>n, (4.1) 
W = a(m/2)i/2, 

has been adopted. By varying the parameters a and /3, the function can be varied 
from monotonically decaying, as in Pataki and Greengard's test [11], to oscillatory. 
Figure 4.1. The r™ term is required for validity of the solution and introducing the 
terms in Tmax scales the amplitude of the cosine on its maximum, so that the maximum 
amplitude of is one, reducing errors caused by very large values of r"*. 

In testing the algorithm, a = 1, P = 16, and R = 16. The parameters varied 
were N the number of blocks in r, M the number of DHT nodes, /? the frequency 
of M^"^ and K. The order n was tested up to 64, with rn = n in the evaluation of 
u'^"^(r). Tests were conducted for direct solution on the DHT grid, to evaluate the 
performance of the underlying method, and with interpolation from the Chebyshev 
nodes. The error measure is the Loo norm: 

max jiti"'' (r, k) — u'"-' (r) | 
max |u(")(r)| 

where ui"'' (r, k) is the computed solution. For clarity, errors greater than 1 have been 
set to 1 on the plots. 

4.1. Solution on DHT nodes. Figure 4.2 shows the error in the solution when 
the input is specified directly on DHT nodes, given as a function of the number of 




(4.2) 
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Fig. 4.2. Error in computation on Hankel transform nodes: solid line /3 = 0; dashed line /3 = 8; 
long dashed line (3 = 16. Left hand column: k = 16; right hand column: k = 1024. FYom top to 
bottom: n = 0, 16, 32, 64. 



nodes M, for /3 = 0,8,16, k ^ 16,1024, and n = 0,16,32,64. In each plot, the 
error behavior is quite similar. For (3 — 0, the solution is not osciUatory and there 
is no aliasing error in the calculation. Thus, the error drops quickly with increasing 
M , as the node density increases, and reaches a minimum, machine precision, around 
M = 64. In contrast, the error for the oscillatory solutions, (3 = 8, 16, remains roughly 
constant for small node number, before dropping quickly to machine precision. The 
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Fig. 4.3. Error in computation on Chebyshev nodes, k = 1024; solid line /3 = 0; dashed line 
/3 = 8; long dashed line /3 = 16. Left hand column: M = 128; right hand column: M = 256. Prom 
top to bottom: n = 0, 16, 32, 64. 



initial failure of the error to reduce with M can be ascribed to aliasing as there are 
insufficient points to capture the oscillatory nature of the input, and the reduction in 
error with M occurs when there is no longer aliasing and the error is controlled by 
the interpolation method. 
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Fig. 4.4. Computation time for solution with n = 64, k = 1024, /3 = 16. Left hand plot uses 
direct calculation at DHT points; circles: computation time; solid line: fit. Right hand plot 
shows computation time for Chebyshev nodes; circles: computation time; curves: linear fit; solid 
line: M = 64; dashed line: M = 128; long dashed line M = 256. 



4.2. Solution on Chebyshev nodes. Figure 4.3 shows similar data, but for 
the error incurred when interpolating from Chebyshev grids onto the DHT nodes. 

The parameters varied arc the same as in Figure 4.2 but error is now plotted as a 
function of the number of nodes NP with M = 128 and 256. Referring to Figure 4.2, 
the choice M = 128 gave machine precision accuracy for ^ = and /3 = 8, though 
not for (3 =- 16, while M = 256 gave machine precision errors for all three values of /3. 

Figure 4.3 shows the result of these choices. The left hand column, where M = 128 
shows the /3 = and /3 = 8 cases reaching the minimum error, though with the error 
controlled by NP. As the node density increases, the interpolation scheme reduces 
the error in the input until the error fixed by the value of M is reached. For /3 = 16, 
the error set by M is greater than machine precision, and the Chebyshev interpolation 
scheme cannot reduce it below the value reached at NP ss 512. Referring to the right 
hand column, all three cases show steady reduction in e down to machine precision, 
but, as might be expected, the /3 = 16 case requires more nodes in order to reduce 
the error to the value set by M. 

4.3. Computation time. Finally, Figure 4.4 shows the computational time for 
the two approaches. The left hand plot shows the time required for direct solution 
on the DHT nodes, as a function of M . Since the DHT is implemented as a matrix 
multiplication, the solution time is expected to scale as and, indeed, the curve 
fit to the data points shows a computation time 0(M^ °°). In solving on a different 
set of nodes, such as the Chebyshev points used here, the number of input nodes 
NP will be greater than M, for reasonable interpolation accuracy. In this case, the 
computation time can be expected to scale as MNP, i.e. proportional to the number 
of nodes and to the number of DHT coefficients used in computing the solution at 
each node. This estimate is borne out by the right hand plot in Figure 4.4 where the 
fits to the computation time have an exponent equal to 1, to three significant figures, 
and the slopes of the lines are proportional to M . 

5. Conclusions. A method for the solution of a modified Bessel equation which 
arises in the solution of Poisson's equation in cylindrical geometries has been pre- 
sented, based on the Hankel transform. Numerical testing has shown the method 
to be accurate over a wide range of wave numbers and orders. We conclude that 
for a proper choice of mesh densities and Hankel transform order, the method can 
achieve machine precision accuracy, for oscillatory and decaying solutions. A sample 
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implementation of the algorithm, written in C, is available from the author. 
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